Different fringe patterns have been shown, and the effect of the wavelength calibration on the spectra was illustrated via the etalon, extended source, and internal calibration source observations. The fringe patterns however are shown "on top" of the spectral baseline of the sources. Different sources have different spectral baselines. Here we look for a way to disentangle the fringe pattern from the source intrinsic spectrum, for the different types of sources considered. In this notebook we look at the etalon data first. The concepts described apply much the same to the broader exercise of calibrating the MRS fringes from extended and point sources.
import funcs
import mrsobs
from numpy import (flipud,mean,diff,array,linspace,pi)
from scipy.optimize import curve_fit
from scipy.interpolate import InterpolatedUnivariateSpline
from matplotlib import pyplot as plt
%matplotlib notebook
We load the images for one band of the MRS:
Additionally the pixel-to-wavelength calibration map and the pixel-to-along-slice position map are imported.
# Define paths to data
workDir = '/Users/ioannisa/Desktop/python/miri_devel/'
cdpDir = workDir+'cdp_data/'
d2cMapDir = workDir+'distortionMaps/'
lvl2path = workDir+'FM_data/LVL2/'
# Get data
band = '1A'
etal_source_sci,etal_source_bkg = mrsobs.FM_MTS_800K_BB_extended_source_through_etalon(lvl2path,band,etalon='ET1A')
# Get wavelength calibration pixel map
d2cMaps = funcs.load_obj('d2cMaps_band{}'.format(band),path=d2cMapDir)
lambdaMap = d2cMaps['lambdaMap']
alphaMap = d2cMaps['alphaMap']
nslices = d2cMaps['nslices']
We subtract background exposures where available
# perform transform
etal_source_bkgsubtr = etal_source_sci-etal_source_bkg
We perform an even-odd row signal correction to the data (caused by the read-out pattern of MIRI detector pixel rows)
etal_source_oddevencorr = funcs.OddEvenRowSignalCorrection(etal_source_bkgsubtr)
What happens when a blackbody spectrum is observed through a Fabry-Pérot resonator (a.k.a. an etalon)? The observation of the MTS extended source seen through etalon ET1A shows an example of what happens. Where there is constructive interference one sees peaks in the spectrum. Where there is destructive interference one sees noise. Important to take note:
In order to compare the etalon spectrum to the transmittance function, the former is normalized by:
# Pixel trace in MRS slice
ypos,xpos = funcs.detpixel_trace(band,d2cMaps,sliceID=nslices/2,alpha_pos=0.)
# Etalon spectrum peaks
etal_source_peaks = funcs.find_peaks(etal_source_oddevencorr[ypos,xpos], thres=0.3, min_dist=6)
# Fit etalon lines
fitparams,errors,fitting_flag,range_ini,range_fin = funcs.fit_etalon_lines(lambdaMap[ypos,xpos],etal_source_oddevencorr[ypos,xpos],etal_source_peaks,fit_func='gauss1d',sigma0=0.002)
lineheights = funcs.get_amplitude(fitparams=fitparams,fitting_flag=fitting_flag)
linecenters = funcs.get_linecenter(fitparams=fitparams,fitting_flag=fitting_flag)
# Take univariate spline through the fitted peaks of the etalon lines
interpolator = InterpolatedUnivariateSpline(linecenters,lineheights,k=3,ext=3)
etal_source_peakprofile = interpolator(lambdaMap[ypos,xpos])
# let's look at the result
plt.figure(figsize=(12,4))
plt.plot(lambdaMap[ypos,xpos],etal_source_oddevencorr[ypos,xpos])
funcs.plot_etalon_fit(fitparams=fitparams,fitting_flag=fitting_flag)
plt.plot(linecenters,lineheights,'ro')
plt.plot(lambdaMap[ypos,xpos],etal_source_peakprofile)
plt.xlabel('Wavelength [micron]',fontsize=20)
plt.ylabel('Normalized signal',fontsize=20)
plt.tick_params(axis='both',labelsize=20)
plt.tight_layout()
plt.figure(figsize=(12,4))
plt.plot(lambdaMap[ypos,xpos],etal_source_oddevencorr[ypos,xpos]/etal_source_peakprofile)
plt.hlines(1,4.88,5.77,linestyle='dashed')
plt.ylim(0)
plt.xlabel('Wavelength [micron]',fontsize=20)
plt.ylabel('Normalized signal',fontsize=20)
plt.tick_params(axis='both',labelsize=20)
plt.tight_layout()
Information about the optical properties of an etalon can be attained by studying the etalon lines in wavenumber space. Assuming that the etalon is plane-parallel, and that the incidence angle is small, subsequent peaks in transmission occur at a constant distance from one another. This distance relates to the optical thickness of the etalon itself.
linecenters_wavelength = linecenters.copy() # microns
linecenters_wavenumber = flipud(10000./linecenters_wavelength) # cm-1
mean_linecenter_separation = mean(diff(linecenters_wavenumber)[1:-1]) # omit first and last data point
# fit straight line through distance data for comparison
popt,pcov = curve_fit(funcs.straight_line,linecenters_wavenumber[1:-2],diff(linecenters_wavenumber)[1:-1])
print r'Mean etalon line separation in wavenumber space is: Δσ = {} cm-1'.format(round(mean_linecenter_separation,2) )
print r'The optical thickness D of the etalon is related to Δσ as: D = 1/(2Δσ) = {} um'.format(round(10000./(2*mean_linecenter_separation),2) )
# Let's look at the results
fig = plt.figure(figsize=(12,5))
axs1 = fig.add_subplot(111)
axs2 = axs1.twiny()
axs1.plot(linecenters_wavenumber[1:-2],diff(linecenters_wavenumber)[1:-1])
axs1.plot(linecenters_wavenumber[1:-2],funcs.straight_line(linecenters_wavenumber[1:-2],*popt))
axs1.hlines(mean_linecenter_separation,linecenters_wavenumber[0],linecenters_wavenumber[-1],linestyle='dashed')
axs1.set_xlabel(r'Wavenumber $\sigma$ [cm-1]',fontsize=20)
axs1.set_ylabel(r'$\Delta \sigma$ [cm-1]',fontsize=20)
axs1.tick_params(axis='both',labelsize=20)
axs2.set_xlim(axs1.get_xlim())
tickmarks = array([1750.,1800.,1850.,1900.,1950.,2000.,2050.])
axs2.set_xticks(tickmarks)
axs2.set_xticklabels(funcs.tick_function(tickmarks))
axs2.set_xlabel('Wavelength [micron]',fontsize=20)
axs2.tick_params(axis='both',labelsize=20)
plt.tight_layout()
Since we work with a normalized etalon transmission, assuming that the etalon's reflecting mirrors are made out of the same material, a single (wavelength-dependent) intensity reflectivity can be determined for the optical element. This can simply be done by looking at the spectral continuum level and relating that to the etalon finesse value, thus reflectivity. The expression of reflectivity as a function of etalon transmission is a second order equation (the equation has two valid solutions, but only one lies between 0 and 1).
# at lower wavelength end
T_e1 = 0.1
F,R_solution1,R_solution2 = funcs.reflectivity_from_continuum(T_e1)
print ' Finesse = {} \n R1 = {} \n R2 = {} (physically valid solution)\n'.format(F,round(R_solution1,2),round(R_solution2,2) )
R_low = R_solution2
# at higher wavelength end
T_e2 = 0.04
F,R_solution1,R_solution2 = funcs.reflectivity_from_continuum(T_e2)
print ' Finesse = {} \n R1 = {} \n R2 = {} (physically valid solution)'.format(F,round(R_solution1,2),round(R_solution2,2) )
R_high = R_solution2
plt.figure(figsize=(12,4))
plt.plot(lambdaMap[ypos,xpos],etal_source_oddevencorr[ypos,xpos]/etal_source_peakprofile)
plt.hlines(1,4.88,5.77,linestyle='dashed')
plt.hlines([T_e1,T_e2],4.88,5.77,'r',linestyle='dashed')
plt.ylim(0)
plt.xlabel('Wavelength [micron]',fontsize=20)
plt.ylabel('Normalized signal',fontsize=20)
plt.tick_params(axis='both',labelsize=20)
plt.tight_layout()
We use the information derived from the data to model the transmittance function of etalon ET1A.
wavelengths = linspace(4.88,5.77,10000) # um
wavenumbers = flipud(10000./wavelengths) # cm-1
R = flipud(linspace(R_low,R_high,10000)) # [-]
D = 1./(2*mean_linecenter_separation) # cm
phi = pi+0.8 # rad
T_e = funcs.FPfunc(wavenumbers,R,D,phi) # [-]
# Let's look at the result
fig,axs = plt.subplots(2,1,figsize=(12,8))
axs[0].plot(lambdaMap[ypos,xpos],etal_source_oddevencorr[ypos,xpos]/etal_source_peakprofile)
axs[0].plot(wavelengths,flipud(T_e))
axs[0].set_ylim(0)
axs[0].set_xlabel('Wavelength [micron]',fontsize=20)
axs[0].set_ylabel('Normalized signal',fontsize=20)
axs[0].tick_params(axis='both',labelsize=20)
axs[1].plot(lambdaMap[ypos,xpos],etal_source_oddevencorr[ypos,xpos]/etal_source_peakprofile)
axs[1].plot(wavelengths,flipud(T_e))
axs[1].set_xlim(5.503,5.692)
axs[1].set_ylim(0)
axs[1].set_title('<Zoomed view>',fontsize=20)
axs[1].set_xlabel('Wavelength [micron]',fontsize=20)
axs[1].set_ylabel('Normalized signal',fontsize=20)
axs[1].tick_params(axis='both',labelsize=20)
plt.tight_layout()
We have shown how the transmittance function of a Fabry-Pérot resonator can be used to model the etalon ET1A data (or any other etalon data) acquired with the MIRI MRS. This was done by normalizing by a spline taken through the etalon transmission peaks. A clear disadvantage of such an approach is that we are not taking into account the beating seen in the data, in fact we are consciously erasing that signature, via the employed normalization. We retrace our steps a bit here to deal with that issue.
# Peaks in previously-derived spline
npixels = 1024; npeaks = 8
etal_source_peakprofile_peaks = funcs.find_peaks(etal_source_peakprofile, thres=0., min_dist=npixels/npeaks/2)
# Take univariate spline through the fitted peaks <== THIS IS THE REAL SPECTRAL BASELINE
interpolator = InterpolatedUnivariateSpline(lambdaMap[ypos,xpos][etal_source_peakprofile_peaks],etal_source_peakprofile[etal_source_peakprofile_peaks],k=3,ext=0)
etal_source_peakprofile_profile = interpolator(lambdaMap[ypos,xpos])
# Let's look at the results
plt.figure(figsize=(12,4))
funcs.plot_etalon_fit(fitparams=fitparams,fitting_flag=fitting_flag)
plt.plot(linecenters,lineheights,'ro')
plt.plot(lambdaMap[ypos,xpos],etal_source_peakprofile)
plt.plot(lambdaMap[ypos,xpos],etal_source_peakprofile_profile)
plt.plot(lambdaMap[ypos,xpos][etal_source_peakprofile_peaks],etal_source_peakprofile[etal_source_peakprofile_peaks],'g*',markersize=20)
plt.xlabel('Wavelength [micron]',fontsize=20)
plt.ylabel('Normalized signal',fontsize=20)
plt.tick_params(axis='both',labelsize=20)
plt.tight_layout()
plt.figure(figsize=(12,4))
plt.plot(lambdaMap[ypos,xpos],etal_source_oddevencorr[ypos,xpos]/etal_source_peakprofile_profile)
plt.plot(linecenters,lineheights/interpolator(linecenters),'ro')
plt.plot(lambdaMap[ypos,xpos],etal_source_peakprofile/etal_source_peakprofile_profile,'b')
plt.plot(lambdaMap[ypos,xpos][etal_source_peakprofile_peaks],etal_source_peakprofile[etal_source_peakprofile_peaks]/interpolator(lambdaMap[ypos,xpos][etal_source_peakprofile_peaks]),'g*',markersize=20)
plt.hlines(1,4.88,5.77,'orange',linestyle='dashed')
plt.hlines(0.846,4.88,5.77,'k',linestyle='dashed')
plt.ylim(0)
plt.xlabel('Wavelength [micron]',fontsize=20)
plt.ylabel('Normalized signal',fontsize=20)
plt.tick_params(axis='both',labelsize=20)
plt.tight_layout()
# compute minima separation (to determine optical thickness of beating component)
beatingminima_wavelength = lambdaMap[ypos,xpos][etal_source_peakprofile_peaks] # microns
beatingminima_wavenumber = flipud(10000./beatingminima_wavelength) # cm-1
mean_beatingminima_separation = mean(diff(beatingminima_wavenumber))
# fit straight line through distance data for comparison
popt,pcov = curve_fit(funcs.straight_line,beatingminima_wavenumber[:-1],diff(beatingminima_wavenumber))
print r'Mean etalon line separation in wavenumber space is: Δσ = {} cm-1'.format(round(mean_beatingminima_separation,2) )
print r'The optical thickness D of the etalon is related to Δσ as: D = 1/(2Δσ) = {} um'.format(round(10000./(2*mean_beatingminima_separation),2) )
print r'The error on the optical thickness D of the etalon is: -{}um +{}um'.format(round(10000./(2*mean_beatingminima_separation) - 10000./(2*diff(beatingminima_wavenumber)[-1]),2 ),round(10000./(2*diff(beatingminima_wavenumber)[0]) - 10000./(2*mean_beatingminima_separation),2) )
print r'The beating fringe is produced by an etalon with an optical thickness: D = {}+{} um = {} um'.format(round(10000./(2*mean_linecenter_separation),2),round(10000./(2*mean_beatingminima_separation),2),round(10000./(2*mean_beatingminima_separation) + 10000./(2*mean_linecenter_separation),2) )
# Let's look at the results
fig = plt.figure(figsize=(12,5))
axs1 = fig.add_subplot(111)
axs2 = axs1.twiny()
axs1.plot(beatingminima_wavenumber[:-1],diff(beatingminima_wavenumber))
axs1.plot(beatingminima_wavenumber[:-1],funcs.straight_line(beatingminima_wavenumber[:-1],*popt))
axs1.hlines(mean_beatingminima_separation,beatingminima_wavenumber[0],beatingminima_wavenumber[-2],linestyle='dashed')
axs1.set_xlabel(r'Wavenumber $\sigma$ [cm-1]',fontsize=20)
axs1.set_ylabel(r'$\Delta \sigma$ [cm-1]',fontsize=20)
axs1.tick_params(axis='both',labelsize=20)
axs2.set_xlim(axs1.get_xlim())
tickmarks = array([1750.,1800.,1850.,1900.,1950.,2000.,2050.])
axs2.set_xticks(tickmarks)
axs2.set_xticklabels(funcs.tick_function(tickmarks))
axs2.set_xlabel('Wavelength [micron]',fontsize=20)
axs2.tick_params(axis='both',labelsize=20)
plt.tight_layout()
# from the depth of the minima caused by the beating we can determine the reflectivity of the etalon producing the beating fringe
T_e1 = 0.846
F,R_solution1,R_solution2 = funcs.reflectivity_from_continuum(T_e1)
print ' Finesse = {} \n R1 = {} \n R2 = {} (physically valid solution)'.format(F,round(R_solution1,2),round(R_solution2,2) )
R_beatingfringe = R_solution2
D_beatingfringe = 1/(2*mean_beatingminima_separation) + 1/(2*mean_linecenter_separation) # cm
phi = pi+0.8 # rad
T_e_beatingfringe = funcs.FPfunc(wavenumbers,R_beatingfringe,D_beatingfringe,phi) # [-]
# Let's look at the results
plt.figure(figsize=(12,4))
plt.plot(lambdaMap[ypos,xpos],etal_source_oddevencorr[ypos,xpos]/etal_source_peakprofile_profile)
plt.plot(linecenters,lineheights/interpolator(linecenters),'ro')
plt.plot(lambdaMap[ypos,xpos],etal_source_peakprofile/etal_source_peakprofile_profile,'b')
plt.plot(wavelengths,flipud(T_e*T_e_beatingfringe))
plt.hlines(1,4.88,5.77,'orange',linestyle='dashed')
plt.hlines(0.846,4.88,5.77,'k',linestyle='dashed')
plt.ylim(0)
plt.xlabel('Wavelength [micron]',fontsize=20)
plt.ylabel('Normalized signal',fontsize=20)
plt.tick_params(axis='both',labelsize=20)
plt.tight_layout()
How do point source etalon data compare to extended source etalon data?